## Different functional forms

library(tidyverse)
library(fixest)
library(broom)
library(ggthemes)
library(interflex)

## Main data
# source("code/01-importdata.R")

#### Bins estimator ####

## Run models
models_bins <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                         scaled_elec_renew, scaled_fit, scaled_quota,
                         scaled_delegates, scaled_gcg,
                         scaled_cf_total_provided, scaled_cf_principal_provided,
                         scaled_cf_total_received, scaled_cf_principal_received) 
                     ~ count_temp_bin_10_lag1 + count_temp_bin_15_lag1 + 
                       count_temp_bin_20_lag1 + count_temp_bin_25_lag1 + 
                       count_temp_bin_30_lag1 + count_temp_bin_35_lag1 + 
                       count_temp_bin_40_lag1 + annual_W_lag1 |
                       iso3c[year] + year, # trends
                       data = analysis_df %>% 
                       mutate(count_temp_bin_20_lag1 = 1))

list_bins <- as.list(models_bins)

# Create empty object to store outputs
tidied_stacked <- matrix(NA, ncol = 6, nrow = 1) %>% 
  set_colnames(c("term", "estimate", "std.error", "statistic", "p.value", "outcome"))

# Loop through models and stack outputs
for (i in 1:11){ 
  tidied <- list_bins[[i]] %>% tidy() 
  tidied$outcome <- list_bins[[i]]$fml %>% gsub("\\~", "", .) %>% .[2] # outcome
  tidied_stacked <- rbind(tidied_stacked, tidied)
}

## Plot bins estimates for two outcomes
tidied_stacked %>% 
  # Clean up
  slice(2:nrow(.)) %>% 
  set_colnames(c("term", "beta", "std_error", "t_stat", "p_value", "outcome")) %>% 
  filter(term!="annual_W_lag1") %>% 
  # Add CIs
  mutate(ci_lower = beta - 1.96*std_error,
         ci_lower_correction = beta - 2.84*std_error,
         ci_upper = beta + 1.96*std_error,
         ci_upper_correction = beta + 2.84*std_error) %>% 
  # Subset the data for two outcomes
  filter(outcome=="scaled_climate_laws" |
           outcome=="scaled_gcg") %>% 
  # Create index for plotting
  mutate(index = str_sub(term, start = 16L, end = -6L),
         index = as.numeric(index)) %>% 
  mutate(index_offset = ifelse(outcome=="scaled_climate_laws", index-0.5, index+0.5)) %>% 
  mutate(Outcome = case_when(outcome=="scaled_climate_laws" ~ "Climate laws", 
                             outcome=="scaled_gcg" ~ "Climate institutions")) %>% 
  mutate(Outcome = factor(Outcome, 
                          levels = c("Climate laws", "Climate institutions"))) %>% 
  # Plot
  ggplot(., aes(index_offset, beta, shape = Outcome)) +
  geom_segment(aes(x = index_offset, xend = index_offset,
                   y = ci_lower_correction, yend = ci_upper_correction),
               size = 1, color = "#79b321") +
  geom_segment(aes(x = index_offset, xend = index_offset,
                   y = ci_lower, yend = ci_upper),
               size = 2, color = "#456613") +
  geom_point(size = 3.5, aes(shape = Outcome)) +
  annotate("point", x = 20.5, y = 0, size = 3.5, shape = 16) +
  annotate("point", x = 19.5, y = 0, size = 3.5, shape = 3) +
  geom_hline(yintercept = 0, linetype = 1) +
  scale_shape_manual(values = c(16,3)) +
  scale_color_manual(values = c("black","black")) +
  scale_y_continuous(breaks = seq(from = -0.015, to = 0.02, by = 0.005)) +
  scale_x_continuous(breaks = seq(10,40,5),
                     labels = c("<10", "10-15", "15-20", "20-25", "25-30", "30-35", ">35")) +
  labs(x="Days in temperature range (Celsius)", 
       y="Regression coefficient \nand confidence intervals") +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        # panel.grid.major.x = element_line(colour = "grey90"),
        # panel.grid.minor.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))
ggsave("text/bins_outcomes_trends.pdf", units = "in", height = 6, width = 8)

## Plot t-statistic for all bin coefficients

set.seed(25) # geom_jitter() used in plot
tidied_stacked %>% 
  # Clean up
  slice(2:nrow(.)) %>% 
  set_colnames(c("term", "beta", "std_error", "t_stat", "p_value", "outcome")) %>% 
  filter(term!="annual_W_lag1") %>% 
  # Add CIs
  mutate(ci_lower = beta - 1.96*std_error,
         ci_lower_correction = beta - 2.84*std_error,
         ci_upper = beta + 1.96*std_error,
         ci_upper_correction = beta + 2.84*std_error) %>% 
  # Create index for plotting
  mutate(index = str_sub(term, start = 16L, end = -6L),
         index = as.numeric(index)) %>% 
  group_by(index) %>% 
  mutate(median_tstat_index = median(t_stat, na.rm=T)) %>% 
  ggplot(., aes(x = index, y = t_stat, group = index)) +
  geom_hline(yintercept = c(-2.84, -1.96, 1.96, 2.84),
             color = c("#79b321", "#456613", "#456613", "#79b321"),
             size = c(2, 1.5, 1.5, 2)) +
  geom_point(aes(x = index, y = median_tstat_index),
             color = "black",
             size = 3,
             shape = 16) +
  geom_jitter(width = 0.5,
              size = 3,
              shape = 1) +
  annotate("point", x = 20, y = 0, size = 3.5, shape = 16) +
  labs(x="Days in temperature range (Celsius)", 
       y="T-statistic") +
  scale_x_continuous(breaks = seq(10,40,5),
                     labels = c("<10", "10-15", "15-20", "20-25", "25-30", "30-35", ">35")) +
  scale_y_continuous(breaks = seq(-4,6,1), 
                     limits = c(-4.1,6.6)) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        # panel.grid.major.x = element_line(colour = "grey90"),
        # panel.grid.minor.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))
ggsave(filename = "text/bins_tstat_trends.pdf", units = "in", height = 6, width = 8)

#### Cumulative shocks, intensification ####

## Two basic results from the canned interflex binning estimator

# For climate laws
analysis_df %>% 
  mutate(ClimateLaws = scaled_climate_laws) %>% 
  interflex(Y = "ClimateLaws",
            D = "hottest_season_temp_lag1", 
            X = "seasonal_heat_shock_cumsum_lag2", 
            Z = "season_W_lag1",
            FE = c("iso3c", "year"),
            estimator = "binning",
            wald = F,
            vcov.type = "cluster",
            cl = "iso3c",
            Dlabel = "Hottest season temperature \n",
            Xlabel = "Cumulative hottest season heat shocks",
            data = ., na.rm=T)
ggsave("text/interflex_climate_laws.pdf", units = "in", height = 6, width = 8)

# For climate institutions
analysis_df %>% 
  mutate(ClimateInstitutions = scaled_gcg) %>% 
  interflex(Y = "ClimateInstitutions",
            D = "hottest_season_temp_lag1", 
            X = "seasonal_heat_shock_cumsum_lag2", 
            Z = "season_W_lag1",
            FE = c("iso3c", "year"),
            estimator = "binning",
            wald = F,
            vcov.type = "cluster",
            cl = "iso3c",
            Dlabel = "Hottest season temperature \n",
            Xlabel = "Cumulative hottest season heat shocks",
            data = ., na.rm=T)
ggsave("text/interflex_gcg.pdf", units = "in", height = 6, width = 8)
                
## For each outcome

m1 <- interflex(Y = "scaled_climate_laws",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m2 <- interflex(Y = "scaled_ecp_zeros",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m3 <- interflex(Y = "scaled_elec_renew",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m4 <- interflex(Y = "fit",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m5 <- interflex(Y = "quota",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m6 <- interflex(Y = "scaled_delegates",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m7 <- interflex(Y = "scaled_gcg",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m9 <- interflex(Y = "scaled_cf_principal_provided",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m8 <- interflex(Y = "scaled_cf_total_provided",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m11 <- interflex(Y = "scaled_cf_principal_received",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)
m10 <- interflex(Y = "scaled_cf_total_received",
                    D = "hottest_season_temp_lag1", 
                    X = "seasonal_heat_shock_cumsum_lag2", 
                    Z = "season_W_lag1",
                    FE = c("iso3c", "year"),
                    estimator = "binning",
                    wald = F,
                    vcov.type = "cluster",
                    cl = "iso3c",
                    data = analysis_df, na.rm=T)

# Plot the L, M, H indicators from interflex
rbind(m1$est.bin[[1]], m2$est.bin[[1]],
      m3$est.bin[[1]], m4$est.bin[[1]],
      m5$est.bin[[1]], m6$est.bin[[1]],
      m7$est.bin[[1]], m8$est.bin[[1]],
      m9$est.bin[[1]], m10$est.bin[[1]], 
      m11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  group_by(Bin) %>% 
  mutate(n = row_number()) %>% 
  group_by(outcome) %>% 
  mutate(m = row_number()) %>% 
  ungroup() %>% 
  mutate(index = ifelse(m == 1, n-0.2,
                        ifelse(m==3, n+0.2, n))) %>% 
  # Plot 
  ggplot(., aes(beta, index)) +
  geom_segment(aes(x = ci_lower_correction, xend = ci_upper_correction,
                   y = index, yend = index),
               size = 1, color = "#79b321") +
  geom_segment(aes(x = ci_lower, xend = ci_upper, 
                   y = index, yend = index), 
               size = 2, color = "#456613") +
  geom_point(size = 3.5, aes(shape = Bin)) +
  scale_shape_manual(breaks=c("Low","Medium","High"),
                     values = c("L","M","H")) + # 16, 3, 4
  scale_color_manual(values = c("black", "black", "black")) +
  geom_vline(xintercept = 0, linetype = 1, colour = "black") +
  scale_y_continuous(trans = "reverse",
                     breaks = c(1:11),
                     labels = unique(labels_vector$labels)) +
  labs(x="Regression coefficient and confidence intervals", 
       y="") +
  geom_vline(xintercept = 0) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=12))
ggsave("text/interflex_bins_coefplot.pdf", units = "in", height = 8, width = 10)

#### Rolling window ####

r1 <- interflex(Y = "scaled_climate_laws",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r2 <- interflex(Y = "scaled_ecp_zeros",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r3 <- interflex(Y = "scaled_elec_renew",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r4 <- interflex(Y = "fit",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r5 <- interflex(Y = "quota",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r6 <- interflex(Y = "scaled_delegates",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r7 <- interflex(Y = "scaled_gcg",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r9 <- interflex(Y = "scaled_cf_principal_provided",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r8 <- interflex(Y = "scaled_cf_total_provided",
                D = "hottest_season_temp_lag1", 
                X = "seasonal_heat_shock_window_lag2", 
                Z = "season_W_lag1",
                FE = c("iso3c", "year"),
                estimator = "binning",
                wald = F,
                vcov.type = "cluster",
                cl = "iso3c",
                data = analysis_df, na.rm=T)
r11 <- interflex(Y = "scaled_cf_principal_received",
                 D = "hottest_season_temp_lag1", 
                 X = "seasonal_heat_shock_window_lag2", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)
r10 <- interflex(Y = "scaled_cf_total_received",
                 D = "hottest_season_temp_lag1", 
                 X = "seasonal_heat_shock_window_lag2", 
                 Z = "season_W_lag1",
                 FE = c("iso3c", "year"),
                 estimator = "binning",
                 wald = F,
                 vcov.type = "cluster",
                 cl = "iso3c",
                 data = analysis_df, na.rm=T)

# Plot the L, M, H indicators from interflex
rbind(r1$est.bin[[1]], r2$est.bin[[1]],
      r3$est.bin[[1]], r4$est.bin[[1]],
      r5$est.bin[[1]], r6$est.bin[[1]],
      r7$est.bin[[1]], r8$est.bin[[1]],
      r9$est.bin[[1]], r10$est.bin[[1]], 
      r11$est.bin[[1]]) %>% 
  as.data.frame() %>% 
  select(beta = coef, se = sd, ci_lower = CI.lower, ci_upper = CI.upper) %>% 
  mutate(ci_lower_correction = beta - 2.84*se,
         ci_upper_correction = beta + 2.84*se) %>% 
  mutate(Bin = rep(c("Low", "Medium", "High"), times = 11)) %>% 
  mutate(Bin = factor(Bin, c("Low", "Medium", "High"))) %>% 
  mutate(outcome = rep(select(analysis_df, starts_with("scaled")) %>% names,
                       each = 3)) %>% 
  group_by(Bin) %>% 
  mutate(n = row_number()) %>% 
  group_by(outcome) %>% 
  mutate(m = row_number()) %>% 
  ungroup() %>% 
  mutate(index = ifelse(m == 1, n-0.2,
                        ifelse(m==3, n+0.2, n))) %>% 
  # Plot 
  ggplot(., aes(beta, index)) +
  geom_segment(aes(x = ci_lower_correction, xend = ci_upper_correction,
                   y = index, yend = index),
               size = 1, color = "#79b321") +
  geom_segment(aes(x = ci_lower, xend = ci_upper, 
                   y = index, yend = index), 
               size = 2, color = "#456613") +
  geom_point(size = 3.5, aes(shape = Bin)) +
  scale_shape_manual(breaks=c("Low","Medium","High"),
                     values = c("L","M","H")) + # 16, 3, 4
  scale_color_manual(values = c("black", "black", "black")) +
  geom_vline(xintercept = 0, linetype = 1, colour = "black") +
  scale_y_continuous(trans = "reverse",
                     labels = unique(labels_vector$labels),
                     breaks = c(1:11)) + 
  labs(x="Regression coefficient and confidence intervals", 
       y="") +
  geom_vline(xintercept = 0) +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=12))
